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Abstract 

We present a novel approach to handle ring artifacts correction in compressed sens¬ 
ing tomographic reconstruction. The correction is part of the reconstruction process, 
which differs from classical sinogram pre-processing and image post-processing tech¬ 
niques. The principle of compressed sensing tomographic reconstruction is presented. 
Then, we show that the ring artifacts correction can be integrated in the reconstruc¬ 
tion problem formalism. We provide numerical results for both simulated and real 
data. This technique is included in the PyHST2 code which is used at the European 
Synchrotron Radiation Facility for tomographic reconstruction. 

1. Introduction 

1.1. Rings artifaets in tomographie reeonstruetion 

During a tomographic acquisition process, some ffaws in the experimental setup 
can lead to unwanted artifacts. One instance is the presence of defective pixels in the 

PREPRINT: Journal of Synchrotron Radiation A Journal of the International Union of Crystallography 



2 


detector, which appear as lines in the sinogram when the defects are independent of the 
projection angle. These spurious lines give raise to rings artifacts in the reconstructed 
object. Even after preprocessing steps like flat-held correction and median Altering, 
these artifacts can remain and are detrimental to the reconstruction quality. Therefore, 
multiple techniques have been developed to tackle this problem. 

1.2. Related work 

Various techniques have been proposed in the literature to reduce or suppress the 
rings artifacts. As reported in (Rashid et al.^ 2012), these techniques can be classifled 
into two groups : sinogram preprocessing and reconstructed images post-processing. 
The preprocessing methods aims at detecting and correcting the spurious lines in 
the sinogram before applying the reconstruction process, thus, rings do not form if 
the method succeeds. On the other hand, post-processing techniques work directly 
on the reconstructed image, trying to extract the concentric circles and Alter them. 
These methods often perform a transformation into polar coordinates to transform 
the concentric circles into straight lines (Prell et al.^ 2009). 

A comprehensive comparison of ring artifact removal methods can be found in 
(Rashid et al.^ 2012). Although these methods certainly provide satisfactory results 
in their limited framework, the authors report that no existing method is really suit¬ 
able for correcting different types of rings, since they always introduce other distor¬ 
tions. A recent work (E. X. Miqueles & Bermdez, 2014) reports a compressed sensing 
approach for rings artifacts reduction using a total variation denoising of the sinogram 
before calling the reconstruction routine. It is a generalization of Titarenko’s Algo¬ 
rithm (Titarenko et al.^ 2010) which consists in a regularization of the sinogram. This 
can also be classifled in the sinogram pre-processing techniques. 
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2. Preamble 


In this section, we introduce the principle and the formalism of compressed sensing 
tomographic reconstruction. This formalism is extended in section]^ for ring artifacts 
correction. 

2.1. Compressed sensing tomographie reeonstruetion 

Computed tomography aims at reconstructing an image x from a set of projections 
y = P{x). Here y denotes the acquired sinogram, x is the slice to be reconstructed 
and P is the projection operator (while its adjoint P* = P^ is the back-projection 
operator). The classical filtered backprojection algorithm enables to reconstruct the 
image, but the number of projections should be of the same order of the number of 
rows in the image to have an acceptable reconstruction according to the Shannon- 
Nyquist sampling theorem. This is often impracticable, and the subsampling leads to 
artifacts in the reconstructed image. 

Compressed sensing techniques exploit a-priori knowledge on the image, like its spar¬ 
sity, in order to bypass this limitation. Instead of computing a closed form solution like 
in the filtered backprojection technique, tomographic reconstruction by compressed 
sensing amounts to an optimization problem 

X = argmin f{y, x) + g{x) (1) 

X 

where f{y,x) is a fidelity term of x with respect to the acquired data y (henceforth 
f{y^x) is denoted /(x)), and g{x) contains a-priori knowledge on the image. In gen¬ 
eral, the regularization term g{x) makes the problem non-smooth, which precludes 
from using usual Gradient-like algorithms Xn+i = f {xn) • 

Advances in convex analysis provide adapted methods, based on proximal split¬ 
ting methods (Combettes & Pesquet, 2009), which are a generalization of projected 
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gradient. One instance is the Iterative Shrinkage-Thresholding Scheme (ISTA). For 
a functional split into a smooth term / and a non-smooth term the first order 
condition at an optimum x reads 


0 e Vf{x) + dg{x) 

0 e f{x) — X + X + dg{x) 
(Id+a^) (x)e(Id-V/) {x) 

x = (Id+%)"^ (Id-V/) (x) 


which suggests the use of a fixed point iterative scheme. The operator (Id -j-dg) ^ is 
called proximal operator : 


-1 f 1 2 

prox;^^ {x) = (Id+A^fif) {x) = argmin | ^ II® “ ®ll2 + 9{x] 
SO that one step of ISTA reads 

®fc+i = Proxg/L [xk - ^V/(a:fc)^ 


( 2 ) 


where L is the Lipschitz constant of the gradient V/ of the fidelity term. Usually, 
the data fidelity term is a L2 norm, so the calculation of the proximity operator 
of / is straightforward. On the other hand, the proximity operator prox;^^ of the 
regularization term does not always have a closed form expression. 


2.2. Total variation regularization 

Depending on the regularization term used in Q, the reconstruction can yield very 
different results. If the regularization term is null, the problems amounts to a least- 
squares reconstruction. This approach tends to blur the edges in the reconstructed 
slice. A commonly used regularization term for image denoising is the Total Variation 
which was introduced in (Rudin et al.^ 1992) as the Rudin-Osher-Fatemi model. This 
prior has the property to preserve the image edges, which is essential especially in 
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tomographic reconstruction where the object in the sample should be distinguishable. 
In a discrete framework, the Total Variation is the LI norm of the gradient : 


TV(a:) = ||V®||i = 5]^V*)(Ww*)(V (3) 

i 

Given an observed sinogram the total variation tomography reconstruction problem 
aims at finding the regularized image x satisfying 

argmin ~ Px\\‘^ + /3TV (x)| (4) 

where [3 weights the regularization with respect to the data fidelity term, and P is 
the projection matrix. If P is equal to the identity matrix, the problem @ amounts 
to the so-called total variation denoising problem : 


argmin |i ||y - a ;||2 +/3TV (a:)| = prox^TV (?/) (5) 

In (Michel et a/., 2011), a dual approach enables to compute the proximity operator 
of the total variation : 

{ z = argmax i — ||Adivz + y ||2 f 

ll-lloo<i ^ ^ (6) 

X = prox;^TV {y) + ^ div z 

Knowing the proximal map of the Total Variation regularization term, the denoising 
problem can be solved by an iterative shrinkage-thresholding scheme. In PyHST2, we 
use an accelerated version known as the Nesterov algorithm (NESTEROV, 2007) or 
FISTA (Beck & Teboulle, 2009). 

However, tomographic reconstruction is not a simple denoising problem, since a 
projection matrix P appears in the problem Q. Constructing a smooth dual of Q 
would entail to invert P^P which is an ill-posed problem. The total variation tomo¬ 
graphic reconstruction is tackled with two nested FISTA loops, each iteration being 
the solution of a denoising subproblem (Beck &; Teboulle, 2009) 


denoise 




(7) 
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Here L is the Lipschitz constant of V f{x) = P^{Px — y) which is calculated with the 
power method. 

2.3. Dictionary learning 

Total Variation regularization performs well for piecewise-constant images since 
edges and uniform regions are reinforced. However, for non piecewise-constant images, 
the cartoon effect might be prejudicial for the reconstruction quality. Thus, another 
regularization technique has to be considered for such images. 

Most natural images have an intrinsic sparsity which can be recovered by an adapted 
transform or by building the best sparsifying basis. The latter technique is called dic¬ 
tionary learning. Given a set of N acquired signals a number of Nc basis vec¬ 
tors (or atoms) are built. Each signal is expressed as a linear combination 
(rCpp,... of the Nc atoms. Dictionary learning is a joint optimization of the 

atoms D and the coefficients Wp under the constraint of sparsity : 

E ll 2 

\\y - Dwp 

P " ( 8 ) 

s.t. Vp, < S 

Here H-Hq denotes the zero norm counting the number of non-zero component of a 
vector. The dictionary D is learned with K-SVD (Mazhar & Gader, 2008). 

In image processing problems like tomographic reconstruction, the image x is divided 
into N square patches of size mxm pixels. Every patch area of index (p) is expressed 
as a linear combination of the atoms : 

patch(p) = '^k,pV>k (9) 

k 

That is, for pixel i of the image x : 

~ ~ '^Pi) ( 10 ) 

k 

where pi denotes the patch containing the pixel i, and Pp. is the center of this patch 
so that i — Pp. belongs to the patch support. 
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To overcome discontinuities in patches borders, the PyHST2 (Mirone et a/., 2014) 
approach allows the patches to overlap. Following formalism used in PyHST2, lp{i) 
denotes the indicator function of patch (p) while lp(i) denotes the indicator func¬ 
tion of the patch cente]0 Every image pixel is then a linear combination of possibly 
overlapping atoms : 

( 11 ) 

P k 

The tomographic reconstruction problem is 


argmin \\y — Px{W)\\^ + /3dl ll'^lli + P • h{W) 

w 

where W — matrix containing the patches coefficients, and 


( 12 ) 


hpv) = E “ E ^k,p^k{i - rp) 


(13) 




is a term promoting patches overlapping. Again, the optimization problem (12) is 
solved with the FISTA algorithm. 


3. Rings correction in compressed sensing reconstruction 

We now present how rings correction can be handled directly in the reconstruction pro¬ 
cess by integrating additional variables in the functional to minimize. This approach is 
independent of the regularization used and can therefore be applied in various frame¬ 
works like total variation and dictionary learning. Sinogram pre-processing techniques 
modify the acquired data to filter the unwanted lines. This filtering often introduces 
new artifacts. On the other hand, image correction techniques can also add new arti¬ 
facts when circular features are detected as artifacts ; and the forward and backward 
Cartesian-polar coordinate transforms lead to a loss of precision even with a bilin¬ 
ear interpolation. When the rings correction is done in the reconstruction process, the 
^ each patch center is a subset of the patch such that the patch centers are not overlapping 
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data is not modified ; it is an initial guess of the solution for the iterative method. The 
solution minimizes the euclidean distance from the reference slice - the backprojected 
sinogram - under regularization constraints. The rings correction is included in both 
the constraints and the functional to minimize. 

In our approach, the rings correction consists in splitting the sinogram into two 
components : the “genuine” sinogram and the spurious straight lines giving rings 
after back-projection. This splitting is done in the reconstruction process, so the two 
components are updated after each iteration. In the end, only the valid sinogram 
component is kept while the rings variables are discarded. We give two examples of 
frameworks using this approach : total variation regularization and dictionary learning 
reconstruction. 

3.1. Rings correction in total variation framework 

Rings correction can be handled by additional variables r in the fidelity term f{y^ x). 
These artifacts come from a spurious line along projection angles in the sinogram. 
Thus, rings variables are stacked in a vector and added to the sinogram for each 
projection. The fidelity term for one projection reads 

f{x,r)^]^\\y-{Px + r)\\l (14) 

here Px and r do not have the same dimensionality ; Px + r means that a vector of 
rings variables is added to each line of the sinogram as illustrated on Figure [3T| The 
functional F(x,r) is 

F{x, r) = /(x, r) + /3 TV (x) + Pr\\r\\^ (15) 

/3 being a parameter weighting the relative importance of spatial regularization, and 
Pr being a penalization parameter for the rings. 

While sinogram preprocessing techniques filter the lines parallel to the projection 
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angle, this approach forces the sinogram to be decomposed as a sinogram Px and 
rings variables r. The sparsity constraint ll'^lli forces the rings variables to have 
only a few not null components, since the LI norm is a good approximation of the 
sparsity-inducing LO norm. 

We emphasize the fact that the sinogram decomposition into a genuine sinogram 
Px and spurious rings r is not a pre-processing technique ; the rings removal is 
intrinsically part of the reconstruction process. At each iteration, the image x and the 
rings variables r are adapted to minimize the energy F{x^r). 



Fig. 1. Principle of the rings separation. 9 is the projection angle and s is the detector 
bin index. The vertical orange and green lines represent spurious lines giving raise 
to ring artifacts. The decomposition Px -\- r forces the ring values to be captured 
in the vector r (independent of the projection angle). In the end, only the part 
without the rings r is returned. 


The minimum of F{x^r) is found with the total variation regularization solver 


presented in 2.2, This iterative algorithm has a step size 7 = 1/L where L is the 
Lipschitz constant of the gradient V/. This constant is an upper bound of the largest 
eigenvalue of the Hessian V^/. Since the gradient is now taken with respect to both 
image and rings variables, the Hessian is 


'^If 

v5,./ 


-pTp pT 


v?/_ 


P Id 


V^,./ = 

its largest eigenvalue is calculated with the power method. 


( 16 ) 
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Once the Lipschitz constant L is obtained, one iteration of the FISTA denoising 
algorithm needs to compute 


1 , 


Vk+i = proXg/L I Vk-i - -Vf{vk- 


= argmin < - 
^ 2 


z - ( Vk-l - -\/f{Vk-l, 


(17) 


+ i/3TV(z^) + i/3,||z^||i| 


Here v denotes the augmented vector v = containing both image and rings 

variables, Zx (resp. Zr) is the part of vector z containing image (resp. rings) variables. 
Since the squared L2 norm is separable, the proximal operator 0 can be written 


Vk^i = argmin < - 


X - ( Xk -1 - -Vf{xk- 1 , 


+ (^^®) + 2 
+ INr-lli r 


Vk-l - -V/(rfc_i; 


(18) 


1 


= prox^TV/L [Xk-l - j'^xfiXk-l 


In short, the proximal operator is separable with respect to the image and rings 
variables. This is convenient because rings and image variables can be updated by 
solving two separate subproblems in a FISTA iteration. 

The first term in (18) is the denoising problem (©• The second term is the proximity 
operator of the LI norm, which has a closed-form expression called soft thresholding 
operator : 


prox;^l|.||^ {u) — max (|i6| — A , 0) • sign {u) element wise 


(19) 
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3.2. Rings correction in Dictionary Learning framework 


Similarly to total variation (15), the proximal operator is separable with respect to 


patch variables W and rings variables r. The problem (12) becomes 


argmin \\y - {Px{W) + r)\\l + /3dl ll'W^lli + /^r lklli + P ' h{W) 

w 


( 20 ) 


On the other hand, in this case, the non-smooth term is only the LI norm which has 


a closed-form proximal operator (19). Thus, the denoising problem is straightforward. 


and only one FISTA loop is required : 

Xk = soft.threshold;3jjj^/i (^Xk - ^Vwf{xk)j 


Vk = soft _threshold;3^/i ( - - Vr/(rfc) 




i+yrTtf 


( 21 ) 


^k+1 ^k T ( ) {p^k ^k—l) 


rk^i 


4-1 

4+1 


{rk - Vk-i) 


3.3. Preconditioning the optimization problem 


In both cases of Total Variation (15) and Dictionary Learning (3.2), the optimiza¬ 
tion process can be sped up by a precondition process. It is well-known that in the 
back-projected slices, the low spatial frequencies are overrepresented with respect to 
high frequencies. In our case, this leads to an ill-conditioned reconstruction problem. 
Therefore, a ramp filter is applied in the Fourier domain to give the proper weight to 
low and high frequencies. 

The filtering is done before back-projection in every iteration, using a discretized 
version of the high pass ramp filter (Murrell, 1996). If there is only one iteration, the 
optimization then reduces to the standard filtered back-projection. This multiplication 
by a ramp filter can be seen as a preconditioner transforming an ellipsoidal (ill-posed) 
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( 22 ) 


(23) 


problem into a spherical problem, thus fastening the rate of convergence. 

The fidelity term / then becomes 

fix,r) = \\C{y-Px-r)\\l 

where C is the preconditioner. The gradient is 

V f_({CPf{Px-y + r)\ 

Here the operator (CP)"^ — is the filtered back-projection ; so the precondi¬ 

tioner is identified with the high-pass filtering process. We notice that the gradient 
with respect to the rings variables should also be filtered. 

On the other hand, adding rings variables in the functional modifies the condition 
number of the problem. In order to prevent the reconstruction problem to be ill- 
posed, we introduced a multiplicative coefficient a doing a variable substitution : 
S ^ Px + r — y becomes S ^ Px + ar — y. Choosing a enables to pre-condition the 
problem with respect to the ring variables, at the expense of an additional parameter. 


4. Results 


We present here some results for both simulated and real data, and compare our 
method with two mainstream techniques of rings correction : sinogram pre-processing 
based on Fourier-Wavelet de-striping (Miinch et al.^ 2009) and image correction using 
polar coordinates transformation (Prell et a/., 2009). 

4-P Simulated data 

We use the standard test image “Lena” containing both smooth components and 
texture details. The tests are divided in increasing levels of difficulty for rings removal 
methods. First, constant lines are added in the sinogram. These lines independent 
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from the projection angle give raise to rings artifacts in the reconstructed slice. Since 
the spurious lines have a constant value, they should be well handled by pre-processing 
techniques. 

Then, lines with variable intensity are added to the sinogram. This makes the cor¬ 
rection more difficult for pre-processing techniques, especially if the lines have sharp 
variations (i.e high frequencies components). 

Finally, ring-shaped features are added in the phantom before adding spurious lines 
in the sinogram. This case is more challenging because correction methods should 
not remove any feature coming from the phantom (they belong to the “true” image), 
while actually removing rings coming from the sinogram (they come from a flaw in 
the experimental setup). 

The reference sinogram pre-processing technique is the wavelet-Fourier filtering 
(Miinch et a/., 2009). This methods first compute the wavelet decomposition at a 
level L of the sinogram. The vertical detail coefficients Vi at level i G [[1,T]] empha¬ 
size the spurious lines that give raise to rings artifacts. In these coefficients, a spurious 
line is nearly constant along the projection angle, thus, it has only low frequencies in 
the Fourier domain. Filtering these few low frequencies in the Fourier domain enable 
to suppress the line after taking the inverse Fourier transform. The filter used is a 
high-pass Gaussian filter whose standard deviation a tunes the bandwidth. Then, the 
sinogram is reconstructed from these filtered wavelet coefficients. The Matlab imple¬ 
mentation of this method is provided in the author’s article.In the tests, a denotes the 
standard deviation of the Gaussian filter and L is the number of levels of the wavelet 
decomposition. 

The image correction technique used here is Rings Correction in Polar Coordinates 
(RCP) (Prell et a/., 2009). It transforms the image into polar coordinates and performs 
a low-pass filtering in the radial direction. The filtered image is then subtracted from 
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the original image, and a threshold is applied to ignore non-artifact structures. The 
result is filtered in the azimuthal direction. After a transformation into Cartesian 
coordinates, the image should only contain rings artifacts ; these are subtracted from 
the original image. A C++ implementation can be found at (Blair, 2014). In the tests, 
the thresholding parameters are set so that all the image pixels can be considered as 
possibly part of an artifact. The important remaining parameters are the maximum 
ring width W, and the maximum angular arc Oq we expect the rings to have. 

Figure shows the results for the firsts test case. The rings are reduced by the 
sinogram pre-processing technique ([^c), but they do not totally disappear. Besides, 
additional artifacts appear after the correction. The parameters of this method - 
wavelet name, decomposition levels and filter bandwidth - depend on the image type 
and the ring width/intensity. In particular, choosing a wavelet with many vanishing 
moments like in (Miinch et al.^ 2009) does not always yield to better results. The RCP 
performs slightly better ^e) ; a strong artifact is added to the right but the result 
is qualitatively better. The Total Variation regularization entirely removes the rings 

(Us)- 

Figure shows the results for the second test case. The sinogram filtering adds 
many spurious rings Qc). The RCP technique removes most of the rings, but a small 
rings details remain. Again, the difference image shows that strong artifacts are added 
Qf) even if they are not obviously perceptible on the reconstructed image. The Total 
Variation regularization Qg) entirely removes the rings artifacts ; the result is quali¬ 
tatively very close to the original phantom. 

A plot of the rings variables during the total variation reconstruction shows that 
the rings are actually captured by the ring vector r. The six peaks representing the 
detected lines in the sinogram are clearly visible in Figure 


lUCr macros version 2.1.6: 2014/10/01 


15 



Fig. 2. Vector of rings variables for second test case (Figure [^b). The horizontal axis 
goes from zero to the number of bins of the detector - that is, in this simulated 
case, 512 for the 512 x 512 test image. 

Figure shows the results for the third test case. Here two features are added to the 
original phantom ([^a) : a black disk and a semi-circular “ring”. These features are 
part of the phantom, they should not be filtered by rings correction techniques. Lines 
added to the sinogram are not constant along the projection angle, and their width 
can be several pixels. This leads to a back-projected image ^h) with large rings. 
The RCP technique ^ e) is more efficient than the sinogram preprocessing ^c). The 
Total Variation rings removal slightly impacts the semi-circular white feature while 
preserving the black disk. 
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Fig. 3. First test case. 

(a) Original phantom, (b) Result of filtered backprojection after adding constant 
lines in the sinogram, (c) Image back-projected after applying the Munch et al 
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Fig. 4. Second test case. 

(a) Original phantom, (b) Result of filtered backprojection after adding lines of vari¬ 
able width and intensity in the sinogram, (c) Image back-projected after applying 
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Fig. 5. Third test case. 

(a) Original phantom, (b) Result of filtered backprojection after adding lines of vari¬ 
able width and intensity in the sinogram, (c) Image back-projected after applying 
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Beside the visual aspect of the corrected image, we use the Peak Signal to Noise 
Ratio (PSNR) as a measure of the correction quality : 

-j N N 

MSEih,h) = 

i=0 j=0 

where MSE{Ii^l 2 ) is the mean square error between the two images Ii and I 2 of 
dimensions N x N and bit depth of 8 bits. Although PSNR gives a score to the overall 
similarity between the corrected image and the original phantom, it is inconsistent 
with the eye perception of quality. For example, RCP performed better than sinogram 
filtering in these tests, but got a lowest PSNR for cases 2 and 3. The structural 
similarity index gives the same kind of results. 

We were surprised that the sinogram pre-processing technique performed poorly on 
these simulated cases. The reason is probably that the intensity of the lines added 
to the sinogram is too large for this method. The technique presented in (Miinch 
et al.^ 2009) only filters the wavelet detail coefficients : if the sinogram line features 
are too large, they are captured by the approximation coefficients rather than the 
detail coefficients. By trying with lines of smallest amplitude, the Fourier-Wavelet 
method actually worked without adding big artifacts in the reconstructed image. For 
our particular experimental cases, however, the sinogram unwanted lines have non- 
negligible amplitudes, making this method ineffective. 

4^2. Experimental data 

We give here some results for reconstructions performed on real data. The samples 
were kindly provided by the ESRF beamline ID 19. 

4 . 2 . 1 . Syntaetie foam The reconstruction technique was used on a syntactic foam 
sample. In this case, the rings are “large” in the extent that the radius difference 
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between the exterior and the interior of the ring is several pixels. This means that the 
spurious lines in the sinogram have several pixel of width along the detector bins axis, 
forming “bands”. However, the intensities of the lines forming a band has too many 
variations to be efficiently filtered by sinogram pre-processing techniques. This case 
is also difficult for slice-correction algorithms which detect circular features, since the 
sample itself have circular features we do not want to remove. 



Fig. 6. Syntactic foam sample tomography acquired at ESRF ID19, with energy of 
19 keV and pixel size of 0.28/im. (a) Filtered back-projection, (b) Correction with 
RCP technique, using the parameters W — 60 and Oq — 60. (c) Reconstruction 
with Dictionary Learning technique in PyHST2, using the parameters = 0.1, 
/3 dl = 0.035 and p = 20. 
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4-2.2. Rhynie chert We apply the reconstruction on a rhynie chert sample. This sit¬ 
uation is almost the opposite of the previous case : the rings artifacts have a small 
intensity in the reconstructed slice, and the sample borders form a nearly circular 
polygonal shape. These border have a huge amplitude with respect to the rest of the 
sample, and the transition between the border and the interior/exterior is very sharp. 
Thus, slice correction techniques would try to remove the borders before any other 
feature in the slice depending on the thresholding parameters. 

We realized that the rings correction was difficult for the total variation reconstruc¬ 
tion : the procedure added rings tangent to one of the slice borders. It turned out 
that the problem was due to the rotation center for (back)projection improperly set, 
leading to accumulating errors in the iterative reconstruction. Indeed, total variation 
and dictionary learning reconstruction require to compute the projection for the func¬ 
tional, and the back-projection for the functional gradient. If the rotation center for 
these operations is not the same that the one used for actually rotate the sample, 
slight errors appear in the (back)projection ; these error accumulate with the number 
of iterations and take the form of circular features (Figure [^c). 

After setting the correct rotation center, we were able to remove the ring artifacts 
(Figure [^d), especially the one near the center of Figure [^a. In this case, the RCP 
technique did quite well (Figure [^b). 
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Fig. 7. Rhynie chert sample tomography acquired at ESRF ID19, with energy of 
17.6 keV and pixel size of 1.52/im. (a) Filtered back-projection with the correct 
rotation axis, (b) Correction with the RCP technique, using the parameters W = 10 
and Oq — 10. (c) Reconstruction with the Total Variation regularization using the 
incorrect rotation axis, (d) Reconstruction with the Total Variation regularization 
using the correct rotation axis. The parameters were (3 = 3’ 10“^ and (3r = 3 -10~^. 


4.3. Execution time and convergence rate 

In this section we measure the execution time required to obtain an acceptable recon¬ 
struction. All the tests are performed on a machine with an Intel xeon cpu e5-1607 v 2 @ s.ooghz 
processor and a GeForce GTX 750 Ti graphic card. 


lUCr macros version 2.1.6: 2014/10/01 




23 



Fig. 8. Execution time as a function of the number of projection. The image used is 
the 512 X 512 test image “Lena” corrupted with the rings presented in the second 
test case. 


We measured that the execution time is the same with rings correction and without 
rings correction, which is coherent since the ring correction is part of the functional. 
Hence, the execution time gives an insight of the time needed by Compressed Sensing 
reconstruction, but it is not the best choice to measure the influence of rings correction 
on the overall reconstruction time. We now focus on the number of iterations required 
to converge to an acceptable solution. The cost of a single iteration depends on the 
number of projections, as it can be guessed with Figure]^ 

To measure the convergence rate, we use the values of the objective function. For 


Total Variation reconstruction, the objective function is given by (15), it includes 
both the fldelity term (Euclidean distance) and the regularization term (LI norm of 
the image gradient). One can notice that the energy is not a monotonic function of 
the number of iterations, which can be seen as a drawback of FISTA with respect to 
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ISTA. 



Fig. 9. Energy as a function of the number of iterations for the Total Variation 
tomographic reconstruction. The energy is normalized by the energy of the last 
iteration in order to have the same scale in the two cases. The image used is the 
512 X 512 test image “Lena” corrupted with the rings presented in the second test 
case, (a) Evolution of energy with 800 proiections. (b) Evolution of energy with 
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With rings correction, the reconstruction process takes more iterations to converge, 
as illustrated in Figurej^ For simulated and real data, it turned out that a satisfactory 
reconstruction can be achieved with less than 1000 iterations without rings correction. 
When the rings correction is activated, it takes about 2000 iteration to correctly 
remove the ring artifacts. Thus, while the rings correction has no additional cost per 
iteration, it takes nevertheless more iterations to converge to an image with removed 
ring artifacts. The “energy transfer” between the fidelity term \\y — {Px + r )||2 and 
the LI norm of the rings ||r||]^ is actually quite slow. 

The convergence rate also depends on the number of projections. Figurej^shows that 
the reconstruction process converges in 500 iterations (400 with rings correction) for 
200 projections, when it takes about 2000 iterations (1500 without rings correction) 
for 800 projections. It is due to the fact that the weight of fidelity term virtually 
increases as the number of projection increases. To counter-weight this, one has to 
increase the penalty term [3 of the regularization, which makes the energy transfer 
between the fidelity term and the ring variables a little faster. 

The reconstruction quality is still excellent with four times less projections - which 
is a crucial interest of Compressed Sensing tomographic reconstruction. Therefore, this 
method makes sense rather for low-dose tomography when there are few projections ; 
and for now it is quite expensive to use it only for rings correction when a reconstructed 
volume is already available. 

4 . 4 * Conclusions 

We presented a new way to correct the rings artifacts that appear in tomographic 
reconstruction. Compressed sensing tomographic reconstruction is a promising tech¬ 
nique which enables to obtain a slice of good quality with fewer projections. Including 
the rings artifacts correction in the iterative reconstruction process has shown to be 
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efficient while requiring no extra pre or post-processing steps. Besides, additional arti¬ 
facts are less likely to appear thanks to the regularization. This method can be adapted 
to any compressed sensing approach, since the only things to do are modifying the 
functional and the iterative correction step accordingly. 

In a further work, we would like to improve the convergence rate of the rings cor¬ 
rection in order to make it more attractive. We also would like to extend this method 
to lines that are not constant along the projection angle in the sinogram, in order to 
cover more general and difficult cases. 
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